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To address the problem of Landau damping in kinetic turbulence, the forcing of the linearized 
Vlasov equation by a stationary random source is considered. It is found that the time-asymptotic 
density response is dominated by resonant particle interactions that are synchronized with the 
source. The energy consumption of this response is calculated, implying an effective damping rate, 
which is the main result of this paper. Evaluating several cases, it is found that the effective damping 
rate can differ from the Landau damping rate in magnitude and also, remarkably, in sign. A limit 
is demonstrated in which the density and current become phase- locked, which causes the effective 
damping to be negligible; this potentially resolves an energy paradox that arises in the application 
of critical balance to a kinetic turbulence cascade. 



INTRODUCTION 

Though linear Landau damping [IH3] is completely un- 
derstood for an isolated Fourier mode, it is not under- 
stood precisely how this phenomenon manifests in set- 
tings where many modes or degrees of freedom are inter- 
acting chaotically by some nonlinear coupling, i.e. in a 
turbulent setting. It is not clear if Landau's long-time 
solution is useful at all in this context. Indeed, this solu- 
tion corresponds to the most weakly damped root among 
a hierarchy of roots, and only emerges after a transient 
period. For a Fourier mode subject to sustained forcing, 
there is no justification for assuming that it will ever be 
found in this asymptotic state, and thus no reason to sup- 
pose that the Landau damping rate will determine how 
quickly energy is removed from the system. 

In this paper we investigate how linear Vlasov dynam- 
ics respond to a random source. We interpret this source 
as a representation of arbitrarily strong turbulent inter- 
actions. This scenario is in contrast to wave turbulence 
of the "fluid" type, i.e. where the state of the system 
is described by a small number of nonlinearly interact- 
ing fields that support linear wave solutions. The ran- 
dom forcing of the corresponding linearized fluid system 
is typically not a very interesting problem to consider, 
since the response can be anticipated as the resonant re- 
sponse of the wave solutions, and the low-dimensional 
intuition of driven oscillators applies. The problem be- 
comes more interesting when one considers exotic linear 
dynamics such as that of non-normal operators [4] , where 
the eigenspectrum by itself is an incomplete description 
of the linear dynamics [5] . In the case of a collisionless 
plasma, the wave solutions are replaced by an infinite hi- 
erarchy of modes, and the intuition from low-dimensional 
forced linear systems no longer applies. Thus, the prob- 
lem of forced linearized Vlasov equation is inherently in- 
teresting, and, apparently, largely unexplored. Consider- 
ing that the linear solutions given by the Landau roots 
are all exponentially damped, a naive guess would be that 
damping in the nonlinear state could be determined by 
an appropriately weighted average of the Landau rates. 
For the forced Vlasov equation, we find that this guess 



completely misses the mark. 

The calculation made here is a logical step in line with 
the statistical formulation of turbulence, which is not 
concerned with individual solutions of the nonlinear dy- 
namics, but rather seeks to understand generic or univer- 
sal properties of the ensemble. In light of the subtle linear 
dynamics of the Vlasov equation, it should not be surpris- 
ing that the conclusions here are physically nontrivial. It 
is important to emphasize that it is a crucial simplifica- 
tion to assume that the dynamics are essentially random; 
but this assumption clearly also limits the applicability of 
the solution. In particular, the general problem of non- 
linear Landau damping, as addressed in the celebrated 
work of Mouhot and Villani [6], is clearly beyond the 
scope of this approach. 

The results of the calculation are summarized as fol- 
lows. We find that under sustained forcing by a sta- 
tionary random source, Landau damping is supplanted 
by another process, namely the resonant response of the 
particles to the source. This interaction causes the dis- 
tribution function to tend toward coherent steady-state 
solutions that consume energy at a rate different from 
the Landau damping rate. Depending on the frequency 
spectrum of the source and other parameters, the effec- 
tive damping rate (defined in terms of this energy con- 
sumption) can be much smaller, comparable, or much 
larger than the Landau damping rate; in some cases the 
energy can even flow inversely, corresponding to a nega- 
tive effective damping rate. This occurs in the absence of 
linearly unstable eigenmodes, and thus may be thought 
of as an example of "generalized instability" [1] in the 
context of a continuum of stable modes. This is a con- 
text in which generalized instability has been expected by 
some authors to play an insignificant role in the genera- 
tion of turbulence, and thus the findings here constitute 
a possible counterexample to such expectations. 

Our findings shed new light on the phenomenon of 
"critical balance" [7] in kinetic turbulence. Originally 
formulated for magnetohydrodynamic (MHD) turbu- 
lence, critical balance can be extended to a kinetic con- 
text [8|. The basic hypothesis of critical balance is 
that the turbulence cascade will proceed anisotropically 
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through scales where the linear mode frequency matches 
an effective nonlinear turnover frequency. For kinetic tur- 
bulence, however, the presence of Landau damping seems 
at odds with the existence of a cascade at all. That is, 
basic considerations of energy balance imply that a crit- 
ically balanced local cascade, subject to scale- by-scale 
damping (at the rate predicted by linear Landau damp- 
ing), must suffer a loss of energy flux as the cascade pro- 
ceeds to smaller scales [5]. Thus, one might expect the 
energy to dissipate rather than cascade. In this paper, 
we resolve this paradox by demonstrating that critically 
balanced fluctuations can be effectively undamped, even 
when the Landau rate is strong. This may explain recent 
numerical simulations that demonstrate the existence of 
critically balanced cascades [TUHT2"] . 

The solution presented in this work may also help 
point the way to advanced Landau-fluid models [T5] . 
which have been the subject of sustained interest for their 
promise to significantly simplify the analysis of kinetic 
systems. The central question in this context is how to 
model the interactions between a hierarchy of fluid mo- 
ments and how to truncate this hierarchy in a manner 
which retains the important physics of the fully kinetic 
system. The forced solutions of this paper constitute 
exact analytic solutions that can be used to test the ca- 
pabilities of Landau-fluid systems. 

FORMULATION OF THE PROBLEM 

Landau damping, as it was first calculated [1 ], de- 
scribes how plasma perturbations damp according to a 
continuous description (the Vlasov equation), in the ab- 
sence of any model of collisional dissipation. In this form, 
it is a basic example of how irreversible behavior can oc- 
cur in a fundamentally reversible linear system. However, 
discretized models including small but non-zero dissipa- 
tion have been shown to yield the same solution [TJJ [T3] , 
but with the added conclusion that the damping must 
be attributed to the specific physical mechanism of colli- 
sions. By contrast, damping in the (collisionless) Vlasov 
description relies on the continuum to provide a limit- 
less repository for structure in velocity space, with suc- 
cessively smaller and smaller scales forming as damping 
proceeds. It is a comforting fact that these formulations 
ultimately yield the same answer. In this paper, we will 
work within the collisionless continuum description, as it 
is elegant and succinct, and avoids non-universal features 
associated with finite collisionality. 

A popular explanation of Landau damping appeals to 
the image of plasma receiving energy from the electric 
field by "surfing" on waves - with a bit more plasma 
traveling slower than the wave than that traveling faster, 
a net energy transferred from the wave to the plasma, 
thus damping the wave. However, though energy balance 
clearly must be satisfied, the general solution attributed 



to Case [5] and Van Kampen [51 [TB] reveals that linear 
Landau damping is actually due to a systematic smear- 
ing or "phase mixing" of the distribution function that 
is formally equivalent to the free evolution of a popula- 
tion of uncoupled harmonic oscillators with a distribu- 
tion of frequencies. (An alternate but equivalent state- 
ment is that the stable linear system can be transformed 
to action-angle variables [13 HE]-) The analogy between 
oscillator populations and plasmas has fueled fruitful in- 
teraction between the fields [T3] • For nonlinearly coupled 
oscillators, one finds both ordered and disordered states, 
with disordered states exhibiting Landau damping and 
ordered states corresponding to phase-locking or "syn- 
chronization" of the oscillators. It is interesting to note 
that we also find both damping and synchronization in 
the present study, establishing an intriguing qualitative 
connection. 

As an example of the kind of system we would like to 
understand, consider the electrostatic slab drift-kinetic 
equation. This equation describes plasma dynamics in 
the presence of a strong uniform magnetic guide field B = 
zB in a five-dimensional phase space (as, y, z, v\\, v±), 
where v || = z • v and v± = |z X v|. Fourier-transforming 
in position space and defining fcy = z • k we have 

^ + ik n f+i(k n -^)ipF = ^e(k,k>(k')/(k-k'), 

(1) 

where ip is the normalized electrostatic potential satisfy- 
ing nof = 2na J dvnv±dv± f, e(k, k') = eo(z x k') • k is 
the nonlinear coupling coefficient, and w*(k, v±, vn) is a 
linear frequency, which depends on scale lengths of the 
background density rig and temperature T . The dis- 
tribution / represents small plasma fluctuations about 
a large background Fq, satisfying 2ir J dv\\dv±Fo = uq. 
Noting that none of the terms proportional to the dis- 
tribution function have any dependence on v± in their 
coefficients, we can integrate over this variable. We 
then take the nonlinear term as a given function, de- 
noted SQc, v»,t). Finally, as there is no explicit refer- 
ence to k x or k y , we suppress these variables and hence- 
forth only refer to a single scalar wavenumber k = hi 
and velocity variable v = U|| and substitute the notation 
2ir J v±dv ± F -> F (v). 

FORMAL SOLUTION 

The subject of what follows is a very simple equation 
for a single Fourier component of the distribution of par- 
ticles, f(k, v, t). The wavenumber k is now just a param- 
eter and the problem takes the form of a continuum of 
oscillators of natural frequency kv that interact via an 
order parameter, the density n. That is, we are left with 
the one dimensional forced Vlasov equation 
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ikvf + ikc 2 s n{t)G{v) = S(v, t), 



(2) 



where c s is a characteristic wave propagation speed in 
the plasma. The function S(v,t) is the source, which 
we take to represent interaction with a large number of 
other Fourier components that compose a bath of turbu- 
lent fluctuations. The third term on the left hand side of 
Eqn. [l]has been rewritten in terms of a general function 
G(v); we note that the uniform-background (<j* = 0) 
single-species case considered by Landau, it is defined 
G(v) — ~d v F /n . The fluctuation density n(t) is de- 
fined 



where fo(u) = f(u,0). The inversion of our transforma- 
tion is achieved by the formula 



f(v,t) 



duf( U ,t)r(v), 



(8) 



where f u (v) is an eigenmode (e.g. a Case-Van Kam- 
pen mode) of the unforced system (Eqn. [2] with S = 0), 
defined f u (v) = \(u)S(u - v) + c 2 s P[G(v)/(u - v)} and 
P denotes the principal value with respect to the point 
u = v. The eigenmodes are normalized such that X(u) = 



iPj^dvG^/iu 



Plugging our solution into 



this formula and using Eqn. [3] we find the following ex- 
pression for the density: 



(t)= / f(v,t)dv. 



(3) 



Eqn. [2] is a reduced kinetic description in the sense that 
there there is only one velocity dimension and a single 
scalar wavenumber (dependence on the full wavevector k 
is concealed in constants and can also be accounted for 
in the definition of G(v)). 

To solve this equation, we perform an invertible trans- 
formation [17j [18] that is implied by the solution of Van 
Kampen [TS]. (Note that the solution may be sought by 
other well-established methods, e.g. via the Vlasov linear 
response function |20| . or by Laplace transform.) Thus, 
Eqn. [2] is equivalent to 



dt 



+ ikuf = S(u, t), 



where the transformed distribution is defined 
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/+(") 1 /-(") 
D+(u) D_(u)' 



(4) 
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and we define D±(u) = 1 ± 2iric 2 .G±(u). The positive- 
and negative-frequency parts of an arbitrary function 
g(v) are defined in terms of the Fourier transform by 



g±(u) = ± 



±00 



dve tv 



dv- 



2tt 



-g(v)- (6) 



Note that g(u) = g+(u) + g-(u) and D + = D*_ with * 
denoting the complex conjugate. Also, note that the ab- 
sence of the interaction term in Eqn. [4] now makes our 
problem that of forced non-interacting harmonic oscilla- 
tors. This equation may be solved in the frame rotating 
at the oscillator frequency —ku by directly integrating 
with respect to time. We find 



f(u,t) = e- ikut f (u)+e 



^—ikut 



dt'e lkut ' 'S(u,t'), (7) 



n{t) 



du e 



-ikut 



/„(«)+ / dt'e ikut ' S(u,t') 



(9) 

The first term is easily recognized as the solution due 
to Van Kampen of the damping of an initially smooth 
distribution in velocity space. This term tends to zero at 
large times t — ¥ 00 since the factor exp(— ikut) becomes 
increasingly oscillatory, i.e. by the Riemann-Lebesgue 
lemma. We are more interested in the time-asymptotic 
behavior and so we will focus on the second term. If 
we consider the harmonically driven case, i.e. S(v,t) = 
exp(—iQt)S(v), this term yields 



nn(t) 



-ittt 



du 



-i(ku — Vt)t 



i(ku — f2) 



S(u). (10) 



Note that for finite time the integrand is non-singular 
at ku = Q. However, this quantity must be carefully 
evaluated as t — > 00. Defining x = ku — f2, the quantity 
in brackets can be written as 



77774(2;) - ixt{x), 



(11) 



where rj t (x) = sin(xt) / (ttx) and Xt( x ) = (1 ~ cos(xi))/x. 
We can then use the identities lim^oo rj t (x) = S(x) and 
lim^oo X t(x) = P(l) and (i/ir)P J dx'f{x')/{x - x') = 
f+{x) — f-(x) to evaluate the time-asymptotic response 
of the density. Thus for t — > 00 the quantity uq of Eqn. [T0| 
is evaluated as 



lim nn(t) = e 

t— >oo 



, m (2n S + (n/k) \ 

y k D + {n/k) J ' 



(12) 



where we have assumed k > for simplicity. This so- 
lution is a purely oscillatory mode induced by resonant 
forcing of particles with u = Q/k. The secular growth of 
the distribution function at this velocity point dominates 
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over all other contributions at large times. For a general 
forcing function S(v,t) = d£l exp(— ifi,t)S(v, O), we 
can simply integrate this response to obtain 



lim n(t) 



2n 

T 



D+(0/fc) 



(13) 



Note that for a stationary random source, this solu- 
tion is undamped because the spectral components of 
S are uncorrelated, i.e. S is not a smooth function of 
O. Introducing the ensemble average (.), stationarity im- 

plies (§ + (u,n)Sl(u',Cl')^ = 5(0-0') A(u,tt',0). From 
Eqn.[l3]we find the ensemble response 



<M 2 > = 



4t: 2 



00 dQ A(0/fc,0/fc,0) 



\D + (n/kW 



(14) 



This is the ensemble response of the density. Note the use 
of the time-asymptotic solution to compute the ensemble 
response. This is valid if the time-scales present in the 
source are much smaller than the time domain of the 
system, so that the system spends most of its time in the 
asymptotic state. 



To grasp the physical meaning of Eqns. 13 and 14 for 
an actual turbulent system, it may be helpful to consider 
the following Gedankenexperiment. Imagine two boxes of 
plasma, plasma A and plasma B. Plasma A constitutes 
a steady turbulent bath of fluctuations and Plasma B is 
initially uniform and quiescent. Now imagine that we 
have a perfect measurement device capable of exactly re- 
solving the features of the turbulent bath of Plasma A, 
and also a perfect source capable of driving the plasma 
B in an arbitrary fashion. With this experimental setup, 
we measure the complete instantaneous Fourier spectrum 
of Plasma A and, choosing a specific Fourier mode, ex- 
actly reconstruct the signal S(k, v, t) that is driving that 
mode. Then, we use this signal to set the source as an 
input to the initially quiescent plasma. After a short pe- 
riod of time, i.e. a time comparable to the inverse of 
the characteristic frequency of S(k,v,t), plasma B will 
exhibit precisely one Fourier mode at finite amplitude. 
By measuring this mode and comparing it with the same 
mode in plasma A, we should find that the density mo- 
ment of both will be given by Eqn |13| with possible cor- 
rections due to low frequency contributions to S(k, v, t); 
other moments such as the current and temperature fluc- 
tuations should also agree with expressions analogous to 
Eqn(l3j Thus, for all practical purposes, the mode in 
Plasma B should be a "clone" of the mode in plasma A. 
This correspondence depends on the fact that the influ- 
ence of the initial condition will be lost to conventional 
Landau damping. The accuracy of this reproduction also 
depends on the assumption that the frequency spectrum 
at fixed k is peaked to some degree about a characteris- 
tic frequency. Fortunately, there is evidence that plasma 



turbulence does exhibit characteristic frequencies (see for 
instance [2] [22]), even in a strongly nonlinear state. 

We can now compute the average energy consumption 
rate of our solution. By energy, we refer in this case to 
the quantity \n\ 2 /2, which only differs from the physical 
(electrostatic) energy by a constant. Energy balance is 
found by multiplying Eqn. [2] by n* , integrating over ve- 
locity, ensemble-averaging and taking the real part. By 
stationarity we have d t (\n\ 2 ) /2 = and what remains is 



Re[ik (jn*)} = Re[ 



dvS(v) 



(15) 



where we have defined the current j = vdvf. This 
equation expresses the average balance of energy input by 
the source and consumption by the linear phase-mixing 
term (i.e. wave-particle interaction). We define the en- 
ergy input £ = Re[^rt* dvS(v)^], which we can eval- 
uate using Eqn. fT3l 



2tt 

T 



dfl / dv Re[ 



n(0/fc,t>,0) 

D+(n/k) * 



(16) 



where we define IT via (^S + (u, Q)S*(v, 0')y = <5(0 — 

0')II(w, v, O). Finally, let us define the effective damping 
rate implied by Eqns. [l4| and |l6} 



7eff 



(H 2 )' 



(17) 



EXAMPLES 

The Landau damping rate 7^ is determined by the ve- 
locity dependence of the background distribution func- 
tion F (v) and the ratio of the plasma velocity c s to the 
velocity of typical particles, e.g. vt. The response given 



by Eqn. 14 retains this dependence via the function 



In fact the "dispersion relation" for the Landau roots is 
obtained from the zeros of the analytic continuation of 
D + into the lower half plane. Our solution is also deter- 
mined by the statistics of the turbulence, via the function 
S. Let us now examine some simple cases to show how 
7eff can differ from 7^. 



Equation 17 is an exact result, only assuming station- 
arity, but depends on unknown statistics of the source. 
Let us now take a very simple source, and then con- 
sider how the the result applies to a more general class 
of sources. The goal is to determine generic properties 
of the solution, starting from specific examples that are 
relatively easy to follow. 

The spectrum of an appropriate turbulent bath should, 
for fixed k, be peaked about a characteristic frequency 
O ~ w NL , as discussed above. Thus, for simplicity 
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FIG. 1. Comparison of Landau damping rate with effective 
damping rate for the case of uniform background. The Lan- 
dau rate 71, is thin-dashed, the effective rate 7 c ff is solid and 
a = 2c 2 /vt 2 — 1. Also plotted in thick-dashed (gray) is 
the density response, normalized to fit the figure. The Lan- 
dau root is strongly damped with uj « (1.45 — 0.60i)fcuT- This 
mode corresponds to an ion acoustic wave with an ion-electron 
temperature ratio of 1. 



let us consider a source with just one frequency, S — 
exp(-w m t)f M S Q , where f M = exp(-v 2 /v T 2 ) /(y/irv T ) 
gives Maxwellian velocity dependence and So is a 
randomly phased complex amplitude. We substitute 
A(0/k,0/k,0) = 6(0 - w NL )|/ M+ K L /fc)| 2 |S | 2 and 
n(fl/fc,u,fi) = 5(0 - uj NL )fM+(u NIj /k)fM{v)\So\ 2 and 
find 



= A Re r D +(^/ k ) 1 

7cff 27r Ke[ / M+ (o; NL A) J 



(18) 



were we note that /m+^nlA) = ^(OA^t^t)) where 
Z is the plasma dispersion function and £ = w nl /(/cvt). 

As an example, let us first consider G(v) = 2vJm/vt 2 - 
This is the classic case considered by Landau, corre- 
sponding to a single species plasma with a spatially uni- 
form background. We evaluate D + (w NL //c) = 1 + a[l + 
C-^(C)] where a = 2c 2 /vt 2 - Then the Landau damping 
rate is computed from the zeros of D + (ui/k) in the lower 
half-plane. We plot the effective damping rate compared 
with the Landau rate in Fig. [I] Imagining the plasma as 
a single effective oscillator with natural frequency fc«T, 
we can interpret the plotted quantity 7cff/(fct"r) as the 
inverse of the quality factor of this oscillator. We see 
that for high kv^ /w NL the quality factor tends to a con- 
stant of about 0.9, whereas for low kv^/oj^ the quality 
factor becomes infinite; this is because the density and 
current become phase- locked in this limit. 

Fig. [l] shows that effective damping of the density re- 
sponse can be quite significant, even exceeding the Lan- 
dau rate, but becomes negligible for kv^/cj^ < 0.5. 
Thus, if a cascade occurs that is not significantly damped, 
it seems reasonable to conclude that critical balance, now 
stated as the condition kv^ /w NL < 1, must be satisfied. 



FIG. 2. Comparison of Landau damping rate with effective 
damping rate for several cases of the ITG system. Parameters 
are k = 1 and a = 2c 2 /vt 2 = 1; all cases are sub-marginal, 
r\ < 1 + y/E. Landau rates jl are dashed lines, effective rates 
7 e ff are solid lines. The plots are in order, with r\ increasing 
from top to bottom. 



Indeed, the numerical simulation of tokamak turbulence 
by Barnes et al. [11] (see Fig. 3) exhibits a spectrum of 
fluctuation energy that is dominated by kv^ below the 
effective nonlinear frequency. 

As a second example, let us now take G(v) = 
2/m/vt[x — k{1 + r\(x 2 — 1/2)}], where x — v/vt, and 
k and r] are parameters. This corresponds to the slab 
ion-temperature-gradient (ITG) mode, given the follow- 
ing definitions. We take r\ — L h /Lt (where L n and 
Lt are gradient scale lengths of the background den- 
sity and temperature respectively) and n = k y p/ (kL n ) 
(where p — vt/(2Q c ), O c is the ion cyclotron frequency 
and k y is the wavenumber perpendicular to the directions 
of the magnetic field and background gradients); Finite 
Larmor radius effects are neglected for simplicity. The 
homogeneous s ystem (Eqn. [2] with S = 0) is stable for 
< 77 < 1 + y/2(l + a)K- 2 a- 2 + 1 (see e.g. [23 ). We 
find D + (uW*0 = ! + «[! + CZ(0 + k{(t?/2 - l)Z(Q ~ 

We plot a few cases in Fig. [2j Here we see a simi- 
lar behavior as the homogeneous case but find that the 
effective damping can actually take on negative values 
as 77 approaches the marginal point of the ITG instabil- 
ity. The difference of this behavior from the behavior of 
the Landau solution is striking, and also points to the 
intriguing possibility that the ITG mode could sustain 
turbulence in the absence of linear instability. 



FINAL REMARKS 

Although we have assumed the special case of 
Maxwellian forcing to arrive at Eqn. [l8j the absence of 
damping at low kv^/u^ is actually more general than 
it might appear. In fact, we may just as well make the 
substitution /m(v) — >• F(v) in our forcing, where F(v) 
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is a purely real function that falls off rapidly in v - e.g. 
a Maxwellian function multiplied by a polynomial in v. 
Then F + would appear in place of /m+ in Eqn. 18 Using 
a well known identity (see e.g. [16] h we may write F + = 
(F + iF*)/2, where F*{v) = ir^P dv'F (v')/(v - v') 
is the Hilbert transform of F. Because F is real, the real 
and imaginary parts of F + are F/2 and F*/2, respec- 
tively. Consequently, at large argument, the real part of 
F + goes rapidly to zero while the imaginary part decays 
algebraically; the latter can be verified by a multipole 
expansion of the Hilbert transform. Thus, the asymp- 
totic behavior of F + is qualitatively the same as that of 
/m+ aud so the effective damping under this forcing will 
rapidly go to zero at sufficiently small kvr / oj nl . 

If, on the other hand, we consider a forcing function 
whose complex phase depends on v, finite damping can 
occur at small /cut/wnl- Physically, positive damping in 
this case would result if a relative phase (of the correct 
sign) is favored between the current and density moments 
of the source. However, there is no asymmetry in the 
fundamental equations to cause a preferred phase differ- 
ence - that is, the linear eigenmodes of the autonomous 
system (S — 0) have no relative phase (this reflects the 
time-reversal symmetry of the equation, which is also 
obeyed by the nonlinear system, Eqn. [I]). For a turbu- 
lent cascade, such a phase difference may indeed be rep- 
resented in the range where energy is injected (indeed, 
unstable ITG modes exhibit a phase difference between 
the density and current moments), but there is no reason, 
a priori, to expect this phase difference in the so-called 
inertial range, where it is commonly observed that sym- 
metries of the dynamical equations are satisfied "statis- 
tically" Furthermore, note that if we assume that 
the bath of fluctuations, from which the source S is com- 
posed (see Eqn. [I]), has no statistically preferred phase 
difference, then our solution for / will likewise present no 
phase difference, and thus a system having this property 
is self-consistent. 
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